Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 1 February 2008 (MN M?eX style file vl.4) 



Non-parametric Reconstruction of Cluster Mass 
Distribution from Strong Lensing: Modelling Abell 370 

Hanadi M. AbdelSalam^tl ^ Prasenjit Saha^§ & Liliya L.R. Williams^^ 

^Department of Physics (Astrophysics), Keble Rd., 0X1 3RH, Oxford 
'^Institute of Astronomy, Madingley Rd., CBS OH A, Cambridge 



ABSTRACT 

We describe a new non-parametric technique for reconstructing the mass distribu- 
tion in galaxy clusters with strong lensing, i.e., from multiple images of background 
galaxies. The observed positions and redshifts of the images are considered as rigid 
constraints and through the lens (ray-trace) equation they provide us with linear con- 
straint equations. These constraints confine the mass distribution to some allowed 
region, which is then found by linear programming. Within this allowed region we 
study in detail the mass distribution with minimum mass-to-light variation; also some 
others, such as the smoothest mass distribution. 

The method is applied to the extensively studied cluster Abell 370, which hosts 
a giant luminous arc and several other multiply imaged background galaxies. Our 
mass maps are constrained by the observed positions and redshifts (spectroscopic or 
model-inferred by previous authors) of the giant arc and multiple image systems. The 
reconstructed maps obtained for v4370 reveal a detailed mass distribution, with sub- 
structure quite different from the light distribution. The method predicts the bimodal 
nature of the cluster and that the projected mass distribution is indeed elongated 
along the axis defined by the two dominant cD galaxies. But the peaks in the mass 
distribution appear to be offset from the centres of the cDs. 

We also present an estimate for the total mass of the central region of the cluster. 
This is in good agreement with previous mass determinations. The total mass of the 
central region is M = 2.0 — 2.7 x lO^^M^g/irQ^, depending on the solution chosen. 
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1 INTRODUCTION 

Gravitational lensing operates on all scales and provides the 
best way to reconstruct mass distribution, without any prior 
hypothesis about the cluster dynamics or mass-to-light ra- 
tio, on large scales from 100 kpc to a few Mpc; i.e., from 
the innermost regions to the far outskirts of clusters. Mod- 
elling of clusters with giant arcs directly confirms that their 
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innermost regions are dominated by dark matter and thus 
plays an important role in probing the distribution of the 
dark matter in rich clusters. 

In the present paper, we describe a new method for 
reconstructing the cluster Abell 370 [hereafter y4370] non- 
parametrically, using the observational constraints provided 
by strong lensing. It is similar to the method described 
in Saha & Williams 1997 for galaxy- lenses, and here we 
develop it for cluster-lenses. To our knowledge, our tech- 
nique is the first of its kind. The only ingredients needed 
for our reconstruction are the positions of the multiple im- 
ages, their redshifts, the luminosity map of the lensing clus- 
ter and its redshift. The positions of the images are taken 
as rigid constraints while the luminosity distribution is a 
loose constraint subordinate to the lensing data. The model 
strikingly predicts the observed parameters associated with 
each image, the ellipticities and orientations. The recon- 
structed mass distribution compares favourably with the 
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ROSAT/HRI X-ray map which is also bimodal, and more- 
over almost all major peaks visible in our reconstructed mass 
map coincide with the X-ray peaks. 

Abell 370 is a very rich cluster of galaxies at redshift 
Zci = 0.375 and its centre appears to be dominated by two 
bright cD galaxies, which are, together with their associated 
dark matter, mostly responsible for the lensing. The two 
cD galaxies are visible on both the optical ground-based 
CCD images and the HST WFC-1 image. It was just over a 
decade ago that A370 was first recognised as a lens (Lynds 
and Petrosian 1986). This was confirmed by the redshift 
measurement of the observed giant blue arc, Zarc = 0.724 
(Mellier et al 1988, Soucail et al 1988). 

The giant blue luminous arc, the numerous arc(lets) 
and the multiple images observed in yl370, distinguish the 
cluster and have made it a target for extensive studies, ob- 
servations and modelling (Soucail et al 1987, Narasimha & 
Chitre 1988, Kovner 1989, Grossman & Narayan 1989 and 
Kneib et al 1993 [hereafter K93]). K93, from their superb 
ground-based CCD image, presented clear evidence that the 
giant arc consists of at least three multiple merging images, 
characterising it as a cusp-arc. They presented a fit for the 
cluster-lens, with the giant arc as three merging images, as- 
suming a bimodal mass distribution calibrated with the gi- 
ant arc's redshift and parameterised by the ellipticities and 
orientations of the two dominant cD galaxies. Such simple 
mass models may ignore substructures in the cluster mass 
map, leading to imprecise inversion for determining the red- 
shifts of the background images (Kneib et al 1996). However, 
their model-inferred redshifts for some of the multiple im- 
ages were subsequently confirmed by spectroscopy (J. Beze- 
court, personal communication). 

The remarkable structure of the giant arc in yl370 and 
its eastern kink (see below) suggests that the source is strad- 
dling a caustic and exhibits a higher order catastrophe than 
a cusp catastrophe. Careful inspection of the HST image 
revealed that there are more than two breaks within the 
arc and tiny elongated bright knots or granules along the 
arc, which indicates immediately that the arc is, in fact, 
five merging segments. The K93 model reproduced only the 
central three parts as multiple images while the two other 
segments emerge as single images of part of the source. Small 
et al 1996, from their HST observations of the giant arc, de- 
tected a possible bulge and faint spiral structures visible on 
the eastern kink. Thus, they claim that the source is a late 
spiral, which is consistent with the spectroscopic identifica- 
tion. Another feature visible on the HST WFC-1 image is 
the radial arc R, which comprises two merging images across 
the inner critical curve. 

Previous models for 71370 assume a predefined mass 
distribution for the cluster, usually based on the parame- 
ters of the two dominant cD galaxies, i.e., model A370 as 
a bimodal cluster. The free parameters used in character- 
ising each clump are core radii, ellipticity and orientation. 
Such a rigid way of modelling may put the resultant mass 
distribution into a corner of the model space allowed by 
observations. Thus, simple mass models based only on the 
observed parameters of the two dominant galaxies may be 
inaccurate even though they reproduce the multiple images 
correctly. 

The structure of this paper is as follows. In section 2, we 
describe the basic lensing equations used in modelling the 



cluster potential and the bending angles. Sections 3 and 4 are 
devoted to the cluster mass reconstruction, and to testing 
the robustness of the method, respectively. Explicit use of 
Fermat's principle to reproduce the exact location of the 
multiple images and to test the code is explained in section 
5, which also discusses individually all the image systems 
used in our reconstruction. A final discussion is given in 
section 6. We take S7o = 1 and A = throughout. 



2 THE METHOD 

The main observables in any lensed system are the image 
positions, relative magnifications for multiple images if any, 
the source redshifts, and the lens redshift. If the images are 
resolved, image ellipticities and orientations are also observ- 
able. In this section we develop a method for reconstruct- 
ing the mass distribution of the lens using constraints pro- 
vided by (a) the positions of multiple images, (b) image 
orientations and (c) image ellipticities, given the lens and 
source redshifts. However, the reconstruction carried out in 
the present paper uses only the multiple-image positions. 



2.1 Pixellated mass distributions 

For a source at unlensed angular position (3 the time delay 
is 

Tie) = '-^(^ie^^f 



j^(e')\n\e-e'\<fe'^, 



(2.1) 



where S(0') is the surface mass density in the lens plane at 
a position 6' . Since 2l and -Dl are fixed for a given cluster, 
we may as well work with a scaled time delay 

t(6») = i(6» - /3)2 - I o{e') In \e-e'\ d'e' (2.2) 

Ds n J 

where 

a(0) = ^, EeHt^-^^^L. (2.3) 

Note that the critical density Ecrit as defined in 02.3) is for 
a source at infinity; hence the factor of Di^s/Ds in ( [2.2[ ). 

We now consider a pixellated mass distribution for the 
cluster, with amn denoting the surface density of the mn-th 
pixel in units of Ecrit . If a is the pixel size, the total mass is 



\ ^ 



Defining 



(2.4) 



(2.5) 



with the integral coveri ng o nly the mn-th pixel, we can write 
the scaled time delay (p.2[) as 



Tie) = ^{0 - /3f - ^mr.^mn{0) 



Ds 



(2.6) 



Defining 
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off 



(2.7) 



and using Fermat's principle, i.e images are formed at the 
extrema of time delay Vqt{0) — 0, leads to the lens equa- 
tion 



— f3 = ^ (TmnarrmiO) ■ 

mn 

Further defining 

7mn(0) 
5mn{d) 

and 

H 

mn 



2 V de^ del 



(2.8) 

(2.9) 
(2.10) 
(2.11) 

(2.12) 



puts the inverse amplification matrix in the form 



l-^> a™„i?™n(6>),(2.13) 



where t^^^ — re-^Oy and so on. 

The quantity tpmn{0) is the contribution of the mn-th 
pixel to the potential at 0; similarly a,„„(0) is the mn-th 
pixel's contribution to the bending angle, and the second 
derivatives are contributions to terms in the inverse ampli- 
fication matrix. 

In this paper we will consider square pixels. We have 
also experimented with pixels that are Gaussian circles with 
overlapping tails. Gaussian pixels avoid discontinuities in 
the mass that square pixels imply, but the difference in the 
final results is very small. This is because the bending angle 
integrates once over the mass distribution and the potential 
integrates twice, which tends to wash out any effects of mass 
discontinuities at pixel boundaries. 

Explicit expressions for i/jmn(0) and its derivatives for 
square pixels are given in Appendix A. 



2.2 Constraint equations and inequalities 

Since 4^mn{0) and its derivatives are known functions, the 
lens equation (2.8) and the inverse amplification (2.13) are 
linear 



m (Tmn and f3, i.e. the unknowns that the recon- 
struction method must infer. This linearity renders multiple- 
image positions, orientations and ellipticities into linear con- 
straints on the pixellated mass distribution. 

First consider mult iple images. For each image, we write 
the lens equation (2.8) at the observed and thus get a 
two-component constraint equation. But each source (3 in- 
troduces two extra numbers to solve for, so multiple image 
systems actually supply 2( (images) — (sources)) constraints. 

Orientations and ellipticities of single but distorted im- 
ages can also provide linear constraints. Suppose we have 
an image at observed elongated along position angle 
and consider the inverse amplification matrix in coordinates 
0'x,Oy, rotated by (j). We will have 



'XX 

Ty'x' 



(2.14) 



having the same form as in equation ( 2.13 ) but with Hmn[0) 
replaced by its rotated version: 



(2.15) 



where c stands for cos2(/) and s stands for sin 20. Since (ji is 

known (bec ause m easured from the observed image), 

in equation (2.14) is linear in Omn- Consider an image that 



appears more elongated than what is expected based on the 
intrinsic ellipticity distribution of galaxies, but does not look 
like an edge-on spiral. With a rough estimate of the galaxy's 
redshift, and assuming that galaxy's intrinsic ellipticity is 
aligned with its lensing-induced elongation, we can infer that 
its magnification along the 6^1 direction is at least k times 
that along the perpendicular direction. In such a case we 
can write 

k\Txi^' \ < \ry'y'\. (2.16) 

This becomes a linear constraint on the amn if we can infer 
the image parity from a rough idea of where the critical 
curves are, and thus remove the absolute value signs. In cases 
where we can confidently assert that (say) the magnification 
along 9^1 is at least k in absolute value, we can write 



- 1/fc < T^/^/ < 1/k; 



(2.17) 



and here the parity doesn't matter. Finally, for very large 
distortions where it is clear that the amplification eigenvec- 
tors are along 9^/ and 9yi we can write 

Tx'y' = 0; (2.18) 
here again parity doesn't matter. 



Notice that the constraint (2.16), usually with an equal 



ity sign, is the type of information that weak lensing observa- 
tions provide. Because weak lensing generally produces mild 
ellipticity changes, many galaxy images have to be averaged 
over to suppress the noise due to intrinsic galaxy ellipticities 
and enhance the lensing signal. The same procedure can b e 
applied in our case using average magnifications in (2.16), 



instead of magnification of individual galaxies. This opens 
up the possibility of combining strong and weak lensing data 
in a mass reconstruction. In this paper we will explore only 
the multiple-image constraints and leave ellipticities and ori- 
entations for future work. 



2.3 Producing a mass map 

The various observational constraints above, combined of 
course with 



nn > 0, 



(2.19) 



confine the pixellated mass distribution to some allowed re- 
gion. Since the constraints are all linear, it is straightforward 
to find this allowed region by linear programming. But since 
the number of pixels will in practice far exceed the number 
of constraints, the allowed region will contain a vast family 
of mass distributions, all consistent with the observations. 
To obtain mass maps, we need to add more information. 

There are several standard ways of adding and justify- 
ing the extra information. For example, we could ask for the 
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smoothest mass distribution consistent with the data, or for 
a maximum entropy distribution. But for this problem, a 
different figure of merit seems more appropriate: since one 
of the aims of cluster mass reconstruction is to test how well 
light traces mass, it is interesting to study the mass distri- 
bution that follows the light as closely as the lensing data 
allow. We can think of this as a 'minimum M/L variation' 
mass map. We have found that this criterion by itself tends 
to produce artifacts on small scales, so we include a term 
that tends to smooth over the mass distribution on short 
scales. More precisely, our mass maps minimise subject 
to the linear lensing constraints, where 

ran 

+ ^0. ^ ^ (o'm+l.n+l + + 
ran 

^ra-\-\^n—\ ran y. (2.20) 



60 



In equation ( 2.2C ) , L„i„ is the light associated with the mn- 
th mass pixel, scaled so that imn = 1. The first term 

in A^ thus tends to minimise mass-to-light variations. The 
second term in A^ is a discrete version of e** J'(V'^cr)^; mini- 
mizing the integrated square of second derivative is one way 
of smoothing. As before, a is the pixel size, so e can be in- 
terpreted as a smoothing scale. 

The numerical problem we now have to solve is to find 
the minimum of a quadratic function of the Gran (i.e., A^) 
subject to linear equality and inequality constraints involv- 
ing the Gran- This type of problem is known as quadratic 
programming. If there is a solution, it is unique and subrou- 
tines for finding it are widely available. We used the NAG 
routine E04NFF. The technique is limited by storage rather 
than speed: for A*' pixels the storage required is ~ 2A'^^. So 
a few thousand pixels is the current limit. 



3 THE LENS RECONSTRUCTION 

The cluster A370 (^d = 0.375, corresponding to 73l = 
6.1 kpcarcsec"^), hosts not only the first detected giant arc 
but also several multiple resolved images and weakly dis- 
torted arclets (Fort et al 1988, K93). In this paper we will 
consider only multiply imaged systems, which are in a field 
of 2' X 2'. We take the image positions from the recent HST 
images of A370 (Small et al 1996) and the redshifts from 



K93; these are tabulated in Table 3.1 and illustrated on Fig. 
|l[ The critical surface mass density for sources at infinity is 
5.04 X 10^° h^f^ MQurcsec-'^ . 

The results we describe in this paper use pixel size 
a = 2.1"; experiments with other pixel sizes indicated lit- 
tle sensitivity to a. 
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Figure 1. Multiple images systems in the cluster A370. 



System 


Images 


z (arcsec) 


j/(arcsec) 


AO 


Giant arc 






{z = 0.724) 


PI 


-13.58 


-18.53 




P2 


-8.05 


-20.27 




P3 


-3.40 


-20.76 




P4 


1.75 


-19.30 




P5 


6.11 


-16.96 


B 








(z = 0.806) 


B2 


2.43 


9.41 




B3 


7.76 


8.92 




B4 


-24.62 


9.22 


C 








{z = 0.810) 


CI 


0.68 


-4.75 




C2 


4.27 


-6.31 


R 


Radial arc 






(z = 1.3) 


Rl 


-0.49 


-4.17 




R2 


-0.79 


-5.36 


Al 








(z ~ 1.2) 


SI 


-0.10 


47.72 




S2 


1.46 


47.43 


A13 








{z = 1.7 ±0.2) 


Ql 


53.06 


21.05 




Q2 


52.38 


22.41 



Table 1. Positions of the various sets of image system as taken 
from the HST images. 



3.1 Observational constraints 



As explained in Section 2.2, each image contributes two lin- 
ear constraint equations but each source adds two new quan- 
tities to solve for, and hence the number of constraints on 
the mass distribution is 2( (images) — (sources)). 

• The giant arc AO at zao = 0.724 (Soucail et al 1987), 



with its eastern kink, was modelled by K93 as multiple merg- 
ing images of a single background source. However, the iden- 
tification of the multiple images is somewhat complex and 
uncertain. For example, in K93 only the central part is mul- 
tiply imaged. From inspection of the HST image it seems 
plausible that the giant arc is a five image-system, and this 
is the interpretation we follow in this paper. However, we 
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Figure 2. Luminosity map of the cluster A370. Filled numbered 
circles mark the positions of galaxies used to produce the map. 
The numbering of the galaxies follows Mellier et al 1988. 



found that a three-image interpretation gave a very similar 
lens reconstruction. This system supplies 8 constraints. 

• The B2-_B3 and BA images at zb = 0.806 (J. Bezecourt 
private communication) , which is close to the model-inferred 
value of 0.865 in Kneib et al (1994 [hereafter K94]), are 
identified as multiple images of a single source, and thus 
they provide 4 constraints. 

• The faint pair C1-C2 at zc = 0.810 (K93, K94) con- 
tributes 2 constraints. 

• The recently identified radial arc R a,t zr = 1.3 ± 0.2, 
(Smail et al 1996) on the IfST images consists of two main 
elongated segments, in which a few bright granules separated 
by sub-arcseconds are visible. Again we face here the uncer- 
tainty in image identification, and we take only the position 
of the two main segments. Hence, the radial arc provides us 
with 2 more constraints. 

Most of the image systems considered above are located 
in the innermost regions of the cluster and thus will con- 
strain mostly the region they lie in. We added a few more 
constraints from images further out in the cluster, but still 
lying within in our field. We included the arclets Al and 
A13: 

• The arclet Al at zai ~ 1.2 (K94) is a two-component 
distorted image, 5'1-S'2, and adds 2 more constraints. 

• The arclet A13 at zais = l-7to.i5 (K94) is also a two- 
component distorted image, Ql — Q2, and provides 2 extra 
constraints. 

The total from all the above image systems is 20 con- 
straints. 



3.2 Pixellated light map 

The pixellated Ught distribution of the cluster-lens ^370 is 
based on Mellier et al 1988. We used the apparent magni- 
tudes of 27 individual galaxies, including the two dominant 



Table 2. Fiducial parameters of lensing galaxies used in our re- 
construction. (Note that some galaxies have centroids outside our 
2' X 2' field, but are included here because their wings extend into 
our field.) The numbering and B magnitudes are taken from Mel- 
lier et al. 1989. The two cD galaxies in this cluster are #20 & 
#35. 

FIDUCIAL PARAMETERS OF A370 MEMBERS 



MSFM No. 


x(arcsec) ± 0.2 


3/(arcsec) ± 0.2 


B 


'7 
/ 


23.8 


63.5 


22.22 


8 


lie 
li.o 


63.3 


21.91 


9 


O A 


62.2 


21.25 


10 


—31.1 


65.1 


20.13 


12 


18.5 


43.1 


20.86 


13 


22.8 


53.1 


21.80 


14 


44.3 


48.0 


21.69 


15 


56.0 


67.8 


19.68 


16 


39.3 


51.1 


22.52 


17 


—41.6 


35.2 


21.39 


lo 


—47.2 


20.9 


20.98 


20 


—2.8 


25.5 


20.40 


21 


-1.3 


19.5 


21.01 


23 


56.6 


15.6 


21.07 


28 


20.7 


-5.0 


21.09 


29 


-20.5 


6.7 


22.82 


30 


-30.0 


1.8 


22.19 


31 


-53.4 


5.7 


21.85 


32 


-55.2 


10.8 


20.44 


34 


-39.9 


-4.0 


21.49 


35 


-5.3 


-11.3 


20.93 


36 


-56.4 


-24.2 


21.90 


37 


-8.2 


-18.9 


21.51 


52 


10.4 


25.9 


21.85 


54 


-32.5 


46.5 


21.64 


56 


-19.9 


-24.9 


22.04 


58 


15.4 


-34.6 


21.10 



cD galaxies [see Table ^ , to obtain a smoothed luminosity 
map. Each of the 27 galaxies is replaced by a Gaussian hght 
profile of dispersion 6g = 8" 



Lmn oc >^ 10 ' exp 



j = l 



{Q rnix ) 

2^1 



(3.1) 



where rrij is the apparent magnitude of the jth galaxy and Qj 
is the location of the j-th galaxy. This smoothed luminosity 
is plotted in Fig. Q. 



3.3 Minimum MjL variation: ML model 

Now we have all the necessary ingredients to reconstruct 
the cluster mass distribution that follows the light as closely 
(in the least-squares sense) as the lensing data allow. We 
also have the option of smoothi ng th e mass distribution at 
various small scales e (see Eq. 

the lensing data. Our reconstructed mass distribution 



while still satisfying 
for 

smoothing parameter e — 3.5" is presented in Fig. ^ Fig- 
ure ^ shows reconstructions with other values of e. All the 
reconstructions satisfy the lensing constraints precisely, but 
we found that e = 3.5" is the lowest value that eliminates 
evident small-scale artifacts, and in the rest of this paper we 
will mainly discuss this case. 
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Figure 3. Our best reconstructed mass distribution of the cluster A370. The cluster is clearly bimodal, with the two most massive clumps 
in the central region close to two cD galaxies (marked with +) in this cluster, but slightly closer together than the visible galaxies. The 
dashed rectangle encloses the region where we consider the mass reconstruction to be most robust. The contour lines are annotated in 
steps of 0.25 X Sdit- In all the following mass maps we adopt the same contour levels. 



As in K93 our mass reconstruction shows bimodality 
in the cluster, and this persists on mass maps for differ- 
ent values of e — see Fig. ^ However, our non-parametric 
reconstruction shows details of this bimodahty not seen in 
previous work. First, the centres of the two dominant mass 
clumps do not coincide with the two dominant cD galaxies, 
but are shghtly closer together, roughly along the line join- 
ing the two cDs. Second, the southern clump is much more 
massive than the northern clump, although the northern cD 
galaxy is brighter. Third, the mass distribution reveals extra 
substructure in the innermost regions of the cluster, between 
the two cDs. Some of these extra substructures are not as- 



sociated with hght, but are required by the lensing data, 
while some are in the vicinity of cluster galaxies which did 
not make it into the Mellier et al 1989 classification scheme, 
and so were not included in our isoluminosity profile, but 
again our model predicts their requirement by the lensing 
data. Thus our non-parametric model indicates that there 
is a substantial amount of dark matter in the innermost re- 
gions of yl370, which does not closely follow the distribution 
of light. 

Comparison of ROSAT/HRI X-ray map (as seen in Fort 
& MeUier 1994, MeUier et al 1994) and our reconstructed 
maps reveals very similar morphologies. Reassuringly, it re- 
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Figure 4. The reconstructed mass distribution using various values for tlie smootliing scale factor t. It is very reassuring that the 
bimodality is preserved for the different values of t. 



veals a coincidence of the two main central clumps that 
are associated with the two giant cD galaxies. Moreover, 
on larger scales the X-ray map shows an elongation to- 
wards the region of arclet Al as well as secondary peaks 
in north-east direction, which are well matched by peaks in 
our reconstructed maps. This close similarity between the 
ROSAT/HRI X-ray map for A370 and our reconstructed 
maps, even on larger scales, supports the results of our lens 
modelling. 



3.4 Maximally flat model: MF model 

To see whether minimizing M / L variation was introducing 
large biases into the reconstruction, we generated another 
reconstruction, which we call the MF model, with 



= const. 



(3.2) 



The resulting mass distribution, presented in Fig. ^ is the 
map with the minimum dispersion in mass that still repro- 
duces the image positions. It is reassuring to see that in this 
case, we still recover the main features of the minimum M/ L 
variation mass map, most importantly, the bimodal nature 
of the cluster. Comparing maps of Figs. I and I we see that 
even some smaller features are reproduced in both maps. 
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Table 3. Previously published mass estimates of the central re- 
gion of the cluster compared with our estimates. Our estimates 
refer to the 45" X 65" region within the dashed rectangle in Fig. ^ 
Previous workers have usually given mass estimates within a crit- 
ical radius, so the regions considered are not precisely equivalent. 



Model 


10i4h-iMQ 


Hammer 1987 


2. 


Soucail et al 1987 


2. 


Narasimha & Chitre 1988 


4.2 


Grossman & Narayan 1989 


0.93 


Mellier et al 1990 


2. 


This work: 




ML 


2.3 ±0.3 


MF 


2.7 


RL 


2.2 ±0.2 



Figure 5. The reconstructed mass distribution using umtorm 
Lmn, which reveals the bimodal nature of the cluster. 



3.5 Uncertainties in redshifts 

Of the various multiple image systems in A370, only ^40 and 
B have spectroscopic redshifts; as noted in 



the redshifts 

of the other systems are model-inferred (K93 & K94) . To test 
whether our reconstructions are sensitive to these redshifts, 
and hence to model dependencies in the K93 and K94 work, 
we generated some more ML maps where the model-inferred 
redshifts are swapped around. The swapping zc zr —> zc 
and Zc zr ^ zaw zc both produced maps that pre- 
served bimodality as well as most of the extra substructures 
in inner regions of the cluster; however we noticed that more 
mciss tends to be put in the vicinity of high redshift images. 
The permutation zc — > zais zr ^ zc yielded no feasible 
solution; our interpretation is that if zai3 becomes too low, 
it is impossible to make the A13 images far enough apart 
without inconsistencies with data on other images. 

We thus conclude that our reconstructions are only 
weakly sensitive to uncertainties in the model-inferred red- 
shifts used. 

More tests of the robustness of our reconstructions are 
described in the next section. 



3.6 Mass estimates 

Our reconstructions provide mass estimates for the whole 
field, but as we saw above, the mass is well constrained only 
in the central region, where most of the multiple-image sys- 
tems lie. Table H lists our mass estimates for the central 
45" X 65" region (enclosed by a dashed rectangle in Fig. ^ 
for the ML, MF and mean RL (see next section) cases, along 
with previous estimates from the literature. From our dif- 
ferent reconstructions, we conclude that 

M45"x65" = (2.3 ± 0.3)^50' X 10'* A/0. (3.3) 

The mass of the whole 2' x 2' field is much less well 
constrained by our reconstructions. The ML model gives 




-60 -40 -20 20 40 
X( arcsec) 



20 20 40 60 
X( arcsec ) 



Mo,, 



: 4.5 X 10' 



Figure 6. Some of the reconstructed mass distributions for the 
RL case (see text); in each case the rotation angle is stated. 



while the MF model gives 

Af2'x2' 6.6 X lO"ft^o^M0. 

That the ML estimate should be lower is understandable 
when we think about what the minimum criterion [Eq. 
( [2.20] )] does: it follows the lensing requirements in the re- 
gions where there are multiple images, elsewhere it tends to 
extrapolate using the given Lmn- 



4 TESTING THE ROBUSTNESS OF THE 
METHOD 

In this section we will estimate the robustness of our method. 
We produce an ensemble of reconstructions, as described 
below, and calculate the ensemble dispersion in mass as a 
function of position in the lens plane. 

For each reconstruction in the ensemble, we derive Lmn 
by rotating the positions of the all the galaxy members by an 



© 0000 RAS, MNRAS 000, 000-000 



Mass distribution of Abell 370 9 




Figure 7. Contours of Ao-mnin the lens plane. The mass distri- 
bution is strongly constrained in the central region of the cluster, 
i.e., region with high number density of images. 



angle in the lens plane. Let us call these RL, rotated light, 
reconstructions. Figure ^ shows some of these. As in the ML 
and MF reconstructions, the main features of the mass map 
also appear in all the versions of the RL maps. Thus we see 
the reconstructed mass maps are quite robust, especially in 
the regions well sampled by the observed images. We will 
now quantify this statement 

From the RL mass distributions, we can readily calcu- 
late pixel-by-pixel ensemble means and dispersions. Let us 
define 



(4.1) 



which is just the ratio of ensemble dispersion to ensemble 
mean. It quantifies the uncertainties in our reconstruction. 
Figure |^ shows contours of AtTmn in the lens plane. It clearly 
shows that the reconstructed mass distribution is indeed 
very well constrained (A(Tm„ < 0.25) in the inner most re- 
gions of the cluster, where almost all the multiple images lie. 
This implies that the image positions indeed strongly con- 
strain the mass distribution enclosed by them. Moreover, 
there are two extra regions far from the centre which show 
some contours revealing constrained mass distributions in 
their vicinity. These are the t wo re gions with the additional 
images, Al and A13 (see table 3.1), which we added to con- 
strain the outskirts of the cluster. 

Since the observed multiple images are not uniformly 
sprinkled over the lens plane, the reconstructed mass distri- 
bution is not equally well constrained over the whole lens 
plane; areas with more images are more constrained. In or- 
der to visualise this effect, in Fig. ^we plot the average value 
of A(T,„„ computed within circles of radius 10" at a num- 
ber of arbitrarily placed points in the lens plane, against the 
number of images inclosed by the circle. Clearly, the robust- 
ness of the reconstructed map in the lens plane increases 
with the number density of images. 



Figure 8. A plot of number of images versus Aamn, the standard 
deviation of the reconstructed mass distribution. The solid line 
connects the average points of Aamn for the corresponding num- 
bers of images, and the dashed line is a fit. Basically the figure 
shows that regions with more images are well constrained. 




-20 - 



-40 - 



-60 L 
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Figure 9. The arrival time contours for set B. Filled circles mark 
the positions of the observed B2 — B3 and BA. The scale shown 
is in arcseconds. 



5 DISCUSSION OF THE VARIOUS IMAGE 
SYSTEMS 

In this section we discuss extensively all the image systems 
used in our modelling and also examine their predicted shape 
parameters. We start with the simplest configuration of im- 
ages. 



© 0000 RAS, MNRAS 000, 000-000 



10 H.M. AbdelSalam, P. Saha & L.L.R. Williams 




(b) 



(c) 



Figure 10. Inferred topological configurations, of the isochronal 
contours that pass through a saddle point, for the various image 
sets considered in this paper. Here (a) represents the three-image 
configuration that corresponds to the B, C Al and Al'i image 
sets, (b) represents the image configuration of the giant arc AO 
and (c) represents the interesting seven image-configuration ob- 
tained for the radial arc system in some reconstructions — in oth- 
ers the lemniscate with two maxima and a saddle point is replaced 
by a single maximum. 
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5.1 B-System 

A good way to verify that our modelling code actually re- 
produces the observed image positions is to use Fermat's 
principle. In Fig. ^ we show a contour plot of the time delay 
surface for the image set B {zb = 0.806). The extrema and 
saddle points of such a surface mark the locations of the 
images produced by the model and as expected they pre- 
cisely coincide with the observed positions of B2, B3 and 
B4. This fit indeed approves of the model. The topography 
of the arrival time surface for the B image set turns out to 
be a lemniscate — see Fig. px| a below. 

The small angular separation between the nearly merg- 
ing pair B2-B3 and the fact that one is almost a mirror 
image of the other suggests that a critical curve must be 
passing between them. We explored critical curves with a 
recursive code that searches for sign changes, and hence ze- 
ros of det |j4~^(0)|. Figure ^ shows the critical curves for 
the B system, with the images B2, B3 and B4 also indi- 
cated. Mapping the critical curves onto the source plane, 
using the lens equation, gives the caustics. Figure ^ shows 
the caustics for the B system, with the inferred position 
of the B-source also marked. The position of the B source 
with respect to the caustics show that it is experiencing a 
beak-to-beak fold catastrophe. 



Figure 11. The critical curve at zb = 0.806, with the observed 
positions of the images B2, 33 and BA marked by asterisks. 



z = 0.806 




-10 

X (arcsec) 



5.2 C-System, Al and A13 systems 

The time delay surface for the C system reproduce precisely 
the observed position of the pair C1-C2 and also predict a 
third image located on the other side of the lens major axis 
(Fig. |l^). The predicted position of the third image is in the 
same location as predicted before by K93 and coincides with 
a faint blob in the HST image. The topological configuration 
of the isochrones of the C system is again a lemniscate. The 
caustics at the predicted C source position also exhibit a 
beak-to-beak fold catastrophe. 

The arclets Al and A13 follow the same pattern of the 
B and C systems and their topological configuration corre- 
spond to lemniscates too (see Fig. nol-(a)). 



Figure 12. The caustics at zb = 0.806, with an asterisk marking 
the inferred position of the source galaxy B. 

The lensed pair usually denoted as D1-D2, when we 
tried including it, gave results similar to those for the Cl- 
C2 pair. However, we have not included this system in the 
present paper because its redshift is still uncertain (Kneib 
1997, private communication). 

5.3 40-System 

Figure |l^ shows the time-delay contours for the giant arc 
AO. The extrema and saddle points precisely fit the observed 
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Figure 13. Time delay contours for the C-system and tlie arclets 
Al and A13. The filled circles marked CI and C2 show the cor- 
responding observed image positions and moreover a third image 
is predicted for the same pair on the other side of the lens major 
axis. The very top arrival time contour is that of the arclet Al 
and the one on the far right belongs to the arclet A13. Though 
this figure does not show such detail, these arclets are inferred to 
be three-image systems too. 
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Figure 15. The critical curve at 240 = 0.724. The observed 
positions of the five segments are marked with asterisks. 
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Figure 14. Time delay contours for the giant arc AO. Filled cir- 
cles mark the observed position of the five segments of AO. 



positions of the five segments that constitute the giant arc. 
The isochrones of this surface correspond to a lemniscate 
embedded in another lemniscate (Fig. px|-(b)). There is no 
counter-arc. 

Figures 0and|l| show respectively the critical curves 
and caustics at zao = 0.724. The lower part of the latter 
figure reveals a spectacular shape for the cusps. Such cusps 
are called butterfiy cuspoids (Berry & Upstill 1980), and 




-15 -10 -5 5 

X (arcsec) 

Figure 16. The caustics at zao = 0.724. The caustic lines in 
the lower region are shown magnified in an inset. The inferred 
source position, marked by a filled circle, indicates clearly that 
the source is exhibiting a catastrophe which is of a higher order 
than a cusp. 



they are a higher-order catastrophe than just a cusp. The 
predicted position of the A galaxy source is marked by a 
filled circle on Fig. ^ which crosses the caustic at least 
twice. These results for the behaviour experienced by the A 
source galaxy indicates that the giant arc ^0 may be the 
result of higher-order catastrophe than just the elementary 
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Figure 17. Contours of arrival time surface for tlie radial-arc. 
The two filled circle, Rl — R2, mark the position of the two main 
segments considered in our mass modelling. 



Figure 18. The critical lines at the redshift of the observed radial 
arc zn = 1.3. Filled circles mark the position of the two segments 
of the radial arc Rl — R2. 



cusp. If SO, it would be the first giant butterfly-cuspoid arc 
observed. 



5.4 _R-Systein 

Deep HST observations of clusters with spectroscopically 
confirmed giant arcs has revealed lots of detailed substruc- 
tures and lensed features which had not been previously de- 
tected on ground-based CCD images. The exceptional lensed 
feature visible on A370 HST WFC-1 image (Small et al 
1996) is the radial arc R, which does not show up even on 
the deep CFHT images of the cluster (K93). Close inspec- 
tion of the radial arc shows some bright blobs within the two 
main elongated segments. In reconstructing the mass distri- 
bution we used only the constraints given by the position of 
the two main segments, which we call R1-R2. 

^From the time-delay surfaces, we infer that the Rl and 
R2 are two nearby images from a five-image or seven-image 
system. Two of the images are located on either side of the 
lens major axis, while the remainder, including R1-R2, are 
aligned together with small angular separation close to the 
southern cD galaxy. Fig. ^ shows the time-delay surface 
corresponding to Fig. H, and here the radial arc consists 
of five aligned images. With Gaussian pixels, the arc was 
formed out of only three aligned images; the small lemniscate 
was replaced by a single maximum, but otherwise the time 
delay surface (and the mass distribution) was very similar. 
The length of the elongated three- or five-image alignment 
is about 5.2" ± 0.2", which exactly fits the observed length 
of the radial arc. 

For this particular case, we recall that all our recon- 
structed mass distributions predict chaotic mass fluctuations 
between the two cD envelopes and detect the presence of 
dark sub-clumps closer to the southern cD galaxy. These 
massive sub-clumps, if taken individually would have their 



corresponding radial critical curves nearly overlapping (Fig. 
|l^), and hence nearly overlapping radial caustic lines. A sin- 
gle extended source lying in the vicinity of these multi-radial 
caustics (Fig. will lead to an observation of such image 
configuration of very close radial arcs. 

It is of particular interest to notice the dramatic change 
of the critical curves and caustics with redshift. (See for ex- 
ample the sequence of Figures ^9|) . Notice the emer- 
gence of the radial critical curves and caustics at high red- 
shifts (Fig. |l|). These suddenly emerged caustics at high 
redshift are called lips caustics (see Schneider et al 1992). 
Since our reconstruction predicts the R source galaxy is 
straddling these newly formed caustics, we argue that the 
radial arc is a "lips" arc. 

Almost all the parameters associated with the various 
images, e.g. amplification, eccentricity, parity and orienta- 
tion, depend more or less on linear combinations of the 
second derivatives of the projected potential and the im- 
ages redshift. Our model predicts these quantities for each 
individual image and are in good agreement with the ob- 
servations. Table 5.4 lists the orientation angle 6 from the 



horizontal in radians and ei and 62, the two eigenvalues of 
the inverse amplification matrix and their ratio gives the ec- 
centricity. The sign of these eigenvalues specifies the parity 
of the image; e.g., the image is minimum L if both ei and 
62 are positive, is a maximum H if both are negative and a 
saddle S if they have opposite signs. 



6 CONCLUSIONS 

In this paper we have developed a non-parametric technique 
for reconstructing the mass distribution of galaxy clusters 
with strong gravitational lensing. We divide the projected 
lens mass into square pixels, and treat each pixel as an inde- 
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Table 4. Shape parameters of various sets of images as predicted by the model. Almost 
all the parameters, i.e., orientations and ellipticities, predicted by the model are close to 
the observed ones. 
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Figure 19. Caustic lines at zn = 1.3, with the filled circle mark- 
ing the predicted position of the R source. 



pendent contributor to the lensing potential. The observed 
positions of multiple images provide linear constraints on 
this pixellated mass distribution. We can then explore the 
mass distributions that are allowed by these constraints; a 
particularly interesting model is the one that follows the 
light as closely as consistency with the lensing data allow. 

We applied the new method to the well known cluster 
Abell 370. The reconstructed mass maps are invariably bi- 
modal. The two largest mass clumps coincide roughly with 
the two cD galaxies, but are closer together than the visi- 



ble cD galaxies. Also, though the northern cD galaxy ap- 
pears to be brighter, our reconstruction reveals that the 
southern mass clump is much more massive. Our reconstruc- 
tions also show various other features that may be identi- 
fied with features on the X-ray map. In the central region, 
where most of the multiple images are, the mass is well con- 
strained; we estimate the mass of a central 45" x 65" field as 
2.3 ± 0.3 X 10"/i~^Mq. The mass of the whole 2' x 2' field 
is estimated to be 4-6x 10"/i"^ A'/q . 

In addition, we have studied in detail the time delay 
surfaces, critical curves and caustics for the various multiple 
image systems in A370. In particular, we argue that the 
giant arc may be a five-image system at a butterfiy cuspoid 
catastrophe, and that the recently discovered radial arc may 
be part of a seven-image system at a lips catastrophe. 



ACKNOWLEDGEMENTS 

We thank Jean-Paul Kneib for providing us with his 1993 
data, and James Binney for critically reading an early ver- 
sion and for fruitful suggestions. Priyamavada Natarajan, 
Inga Schmoldt and Radek Stompor are gratefully acknowl- 
edged for stimulating discussions. 

HMA acknowledges financial support from the Over- 
seas Research Scheme (ORS) and Oxford Overseas Bursary 
(OOB). LLRW would like to acknowledge the support of the 
PPARC fellowship at the loA, Cambridge, UK. 



REFERENCES 

Berry, M.V. and Upstill, C, 1980. Progress in Optics, 18, 257, 
Bezecourt, J., 1997, private communication 
Fort, B., Prieur, J.L., Mathez, G., Mellier, Y., Soucail, G 
A. & A., 200, L17. 



© 0000 RAS, MNRAS 000, 000-000 



14 KM. AbdelSalam, P. Saha & L.L.R. Williams 



Fort, B. & Mellier, Y., 1994. A. & A. Rev., Volume 5, NO. 4, 
637. 

Grossman, S.A. and Narayan, R., 1989. APJ, 344, 637. 

Hammer, F., 1987. in "Third lAP Astrophys. Meeting on 
High Redshift Galaxies", J. Bergeron, D, Kunth, B. Rocca- 
Volmerange and J. Tran Tranh Van{eds.). Frontieres 1988. 

Kneib, J.-P., Mellier, Y., Fort, B. and Mathez, G., 1993. A. & 
A., 273, 367. 

Kneib, J.-P., Mathez, G., Fort, B., Mellier, Y., Soucail, G. and 

Longaretti, P.-Y., 1994. A. & A., 286, 701. 
Kneib, J.-P., Ellis, R.S., Smail, I.R., Coucli, W., Sharpies, R., 

1996. APJ, 471, 643. 
Kovner,I.,1989. APJ, 337, 621. 

Lynds, R. and Petrosian, V. 1986. Bull. Am. Astronomical Soc., 
18, 1014. 

Lynds, R. and Petrosian, V. 1989. APJ, 336, 1. 

Mellier, Y., Fort,B., Bonnet, H., Kneib, J.-P., 1994, in "Cosmo- 
logical Aspects of X-ray clusters of galaxies" NATO ASI series 
C 441, Seitter W.C. (eds), kluwers, Dordrecht. 

Mellier, Y., Soucail, G., Fort, B. and Mathez, G. 1988. A. & A., 
199, 13. 

Mellier, Y., Soucail, G., Fort, B. , Le Borgne J.-F., Pello, R., 

1990. in Gravitational Lensing, Mellier, Y., Fort, B., Soucail, 

G. (eds), Berlin, Springer. 
Narasimha, D. & Chitre, S.M. 1988. APJ, 322 , 75. 
Saha, P. & Williams, L.L.R., 1997. MNRAS, 292, 148. 
Schneider, P., Ehlers, J. & Falco, E.E. 1992. Gravitational lenses, 

Berlin, Springer- Verlag. 
Smail, I., Dressier, A., Kneib, J-P., Ellis, R.S., Couch, W., 

Sharpies, R. and Oemler, A.Jr., 1996. APJ, 469, 508. 
Soucail, G., Fort, B., Mellier, Y. and Picat, J.-P., 1987. A. & A., 

172, L14. 

Soucail, G., Mellier, Y., Fort, B. Mathez, G. and Cailloux, M., 
1988. A. & A., 191, L19 



d 



^mn(a:,2/) 



1 



y arctan — ^ ln{x^ + y'^) ~ x \ , 



dy TT y y 

and the second derivatives via the indefinite integrals 
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APPENDIX A: 

This Appendix gives expre ssio ns for the integral over indi- 



vidual pixels in equation (2.5), and its derivatives. In the 



following we write x, y for 6x,6y. 

The integrals over pixels are most conveniently evalu- 
ated by first considering the corresponding indefinite inte- 
grals and then differencing between pixel-corner values. The 



indefinite integral correspond to equation (2.5) is 



7 , , 1/2 y 2 X 

•ipmn(x,y) — -TT-XX atctau — hy arctan — 
ZTT \ x y 



-\-xy\n{x +y 



Taking the pixel size as a and noting that the mn-th pixel 
is centred at {ma,na), the desired definite integral is 

V'mn(a;,j/) = 'ilimn {x - (m + ^)a,y - {n + ^)aj 

+ i)mn {x ~ {m ~ \)a,y - [n - ^)a) 

- i!mn.{x ~ (m^ \)a,y - {n - \)a) 

- i>mn {x - {m- i)a, y - (n + i)a) . 

Similarly, we can compute the first derivatives of i^mn{x,y) 
via the indefinite integrals 
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dx 
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